

%% "true" DGP point estimates
load('../true_DGP/tDGP.mat','DGP')
true_beta=DGP(1:(2*N));
clear DGP


%% reset rng seed
rng('default')


%%
NSP=1e+3; % number of simulated sample paths
T_spg=1;
epsilon=mvnrnd(zeros(N,1),Sigma,(T-(23-T_spg-1))*NSP)'; % growth shocks
epsilon=reshape(epsilon,N,T-(23-T_spg-1),NSP);


%% small shock to UK, France, Germany, (1973-T_spg):1973
P_spg=zeros(size(P(:,:,23-T_spg+1:T)));
parfor nsp=1:NSP
    tic
    ep=epsilon(:,:,nsp);
    y_spg=y;
    y_spg([24,28,32],23-T_spg:23)=-0.02+ep([24,28,32],1:T_spg+1);
    PP_spg=zeros(size(P(:,:,23-T_spg+1:T))); %#ok<*PFBNS>
    BB_spg=zeros(2*N,T-23+T_spg);
    % 1973-T_spg
    nm=find(M(:,23-T_spg)==1);
    DD=[diag(1-D(nm,23-T_spg)),diag(D(nm,23-T_spg))];
    PP_spg(:,:,1)=P(:,:,23-T_spg);
    PP_spg([nm;N+nm],[nm;N+nm],1)=PP_spg([nm;N+nm],[nm;N+nm],1)+DD'*(Sigma(nm,nm)\DD);
    Pd_spg=sqrt(diag(PP_spg(:,:,1)\eye(2*N)));
    BB_spg(:,1)=B(:,23-T_spg);
    BB_spg([nm;N+nm],1)=BB_spg([nm;N+nm],1)+PP_spg([nm;N+nm],[nm;N+nm],1)\(DD'*(Sigma(nm,nm)\(y_spg(nm,23-T_spg)-DD*BB_spg([nm;N+nm],1))));
    U0_spg=exp(alpha([30:31,45],1)+theta(1)*(Pd_spg([30:31,45],1)*sgi_nodes(:,1)'+BB_spg([30:31,45],1)+...
        ep_sd([30:31,45],1)*sgi_nodes(:,2)'));
    U0_spg=(U0_spg./(1+U0_spg))*sgi_weights;
    U1_spg=exp(alpha([30:31,45],1)+theta(2)*(Pd_spg(N+[30:31,45],1)*sgi_nodes(:,1)'+BB_spg(N+[30:31,45],1)+...
            ep_sd([30:31,45],1)*sgi_nodes(:,2)')-f([30:31,45],1));
    U1_spg=(U1_spg./(1+U1_spg))*sgi_weights;
    DD_spg=D;
    DD_spg([30:31,45],23-T_spg+1)=M([30:31,45],23-T_spg+1).*(U1_spg>=U0_spg);
    for t=1:T-23+T_spg-1
        if DD_spg(30,23-T_spg+t)~=D(30,23-T_spg+t)
            y_spg(30,23-T_spg+t)=true_beta(30)*(1-DD_spg(30,23-T_spg+t))+true_beta(N+30)*DD_spg(30,23-T_spg+t)+ep(30,t+1);
        end
        if DD_spg(31,23-T_spg+t)~=D(31,23-T_spg+t)
            y_spg(31,23-T_spg+t)=true_beta(31)*(1-DD_spg(31,23-T_spg+t))+true_beta(N+31)*DD_spg(31,23-T_spg+t)+ep(31,t+1);
        end
        if DD_spg(45,23-T_spg+t)~=D(45,23-T_spg+t)
            y_spg(45,23-T_spg+t)=true_beta(45)*(1-DD_spg(45,23-T_spg+t))+true_beta(N+45)*DD_spg(45,23-T_spg+t)+ep(45,t+1);
        end
        BB_spg(:,t+1)=BB_spg(:,t);
        PP_spg(:,:,t+1)=PP_spg(:,:,t);
        nm=find(M(:,23-T_spg+t)==1);
        DD=[diag(1-DD_spg(nm,23-T_spg+t)),diag(DD_spg(nm,23-T_spg+t))];
        PP_spg([nm;N+nm],[nm;N+nm],t+1)=PP_spg([nm;N+nm],[nm;N+nm],t+1)+DD'*(Sigma(nm,nm)\DD);
        Pd_spg=sqrt(diag(PP_spg(:,:,t+1)\eye(2*N)));
        BB_spg([nm;N+nm],t+1)=BB_spg([nm;N+nm],t+1)+PP_spg([nm;N+nm],[nm;N+nm],t+1)\(DD'*(Sigma(nm,nm)\(y_spg(nm,23-T_spg+t)-DD*BB_spg([nm;N+nm],t))));
        U0_spg=exp(alpha([30:31,45],1)+theta(1)*(Pd_spg([30:31,45],1)*sgi_nodes(:,1)'+BB_spg([30:31,45],t+1)+...
            ep_sd([30:31,45],1)*sgi_nodes(:,2)'));
        U0_spg=(U0_spg./(1+U0_spg))*sgi_weights;
        U1_spg=exp(alpha([30:31,45],1)+theta(2)*(Pd_spg(N+[30:31,45],1)*sgi_nodes(:,1)'+BB_spg(N+[30:31,45],t+1)+...
                ep_sd([30:31,45],1)*sgi_nodes(:,2)')-f([30:31,45],1));
        U1_spg=(U1_spg./(1+U1_spg))*sgi_weights;
        DD_spg([30:31,45],23-T_spg+t+1)=M([30:31,45],23-T_spg+t+1).*(U1_spg>=U0_spg);
    end
    D_spg(:,:,nsp)=DD_spg;
    B_spg(:,:,nsp)=BB_spg;
    P_spg=P_spg+PP_spg/NSP;
end
D_spg_ss=1*(mean(D_spg,3)>0.5);
B_spg_ss=mean(B_spg,3);


%% for Figure 8

% mean beliefs
writetable(table(B(N+[30:31,45],:)-B([30:31,45],:)),'../../Figures/Figure_8/spg_original.csv')
writetable(table(B_spg_ss(N+[30:31,45],:)-B_spg_ss([30:31,45],:)),'../../Figures/Figure_8/spg_cfactual.csv')

% belief uncertainty
writetable(table(uncertainty([30:31,45],:)),'../../Figures/Figure_8/spg_unc_original.csv')
Pinv_spg=P_spg;
for t=1:size(P_spg,3)
    Pinv_spg(:,:,t)=P_spg(:,:,t)\eye(2*N);
end
unc_spg=zeros(N,size(P_spg,3));
for n=1:N
    for t=1:size(P_spg,3)
        unc_spg(n,t)=sqrt(Pinv_spg(N+n,N+n,t)+Pinv_spg(n,n,t)-2*Pinv_spg(N+n,n,t));
    end
end
writetable(table(unc_spg([30:31,45],:)),'../../Figures/Figure_8/spg_unc_cfactual.csv')



